{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 第6章Logistic回归与最大熵模型-习题\n",
    "\n",
    "### 习题6.1\n",
    "&emsp;&emsp;确认Logistic分布属于指数分布族。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**  \n",
    "**第1步：**  \n",
    "首先给出指数分布族的定义：  \n",
    "对于随机变量$x$，在给定参数$\\eta$下，其概率分别满足如下形式：$$p(x|\\eta)=h(x)g(\\eta)\\exp(\\eta^Tu(x))$$我们称之为**指数分布族**。  \n",
    "其中：  \n",
    "$x$：可以是标量或者向量，可以是离散值也可以是连续值  \n",
    "$\\eta$：自然参数  \n",
    "$g(\\eta)$：归一化系数  \n",
    "$h(x),u(x)$：$x$的某个函数  \n",
    "\n",
    "----\n",
    "\n",
    "**第2步：**证明伯努利分布属于指数分布族  \n",
    "伯努利分布：$\\varphi$是$y=1$的概率，即$P(Y=1)=\\varphi$  \n",
    "$\\begin{aligned}\n",
    "P(y|\\varphi) \n",
    "&= \\varphi^y (1-\\varphi)^{(1-y)} \\\\\n",
    "&= (1-\\varphi) \\varphi^y (1-\\varphi)^{(-y)} \\\\\n",
    "&= (1-\\varphi) (\\frac{\\varphi}{1-\\varphi})^y \\\\\n",
    "&= (1-\\varphi) \\exp\\left(y \\ln \\frac{\\varphi}{1-\\varphi} \\right) \\\\\n",
    "&= \\frac{1}{1+e^\\eta} \\exp (\\eta y)\n",
    "\\end{aligned}$  \n",
    "其中，$\\displaystyle \\eta=\\ln \\frac{\\varphi}{1-\\varphi} \\Leftrightarrow \\varphi = \\frac{1}{1+e^{-\\eta}}$  \n",
    "将$y$替换成$x$，可得$\\displaystyle P(x|\\eta) = \\frac{1}{1+e^\\eta} \\exp (\\eta x)$\n",
    "对比可知，伯努利分布属于指数分布族，其中$\\displaystyle h(x) = 1, g(\\eta)= \\frac{1}{1+e^\\eta}, u(x)=x$  \n",
    "\n",
    "----\n",
    "\n",
    "**第3步：**  \n",
    "广义线性模型（GLM）必须满足三个假设：\n",
    "1. $y | x;\\theta \\sim ExponentialFamily(\\eta)$，即假设预测变量$y$在给定$x$，以$\\theta$为参数的条件概率下，属于以$\\eta$作为自然参数的指数分布族；  \n",
    "2. 给定$x$，求解出以$x$为条件的$T(y)$的期望$E[T(y)|x]$，即算法输出为$h(x)=E[T(y)|x]$  \n",
    "3. 满足$\\eta=\\theta^T x$，即自然参数和输入特征向量$x$之间线性相关，关系由$\\theta$决定，仅当$\\eta$是实数时才有意义，若$\\eta$是一个向量，则$\\eta_i=\\theta_i^T x$\n",
    "\n",
    "----\n",
    "\n",
    "**第4步：**推导伯努利分布的GLM  \n",
    "已知伯努利分布属于指数分布族，对给定的$x,\\eta$，求解期望：$$\\begin{aligned}\n",
    "h_{\\theta}(x) \n",
    "&= E[y|x;\\theta] \\\\\n",
    "&= 1 \\cdot p(y=1)+ 0 \\cdot p(y=0) \\\\\n",
    "&= \\varphi \\\\\n",
    "&= \\frac{1}{1+e^{-\\eta}} \\\\\n",
    "&= \\frac{1}{1+e^{-\\theta^T x}}\n",
    "\\end{aligned}$$可得到Logistic回归算法，故Logistic分布属于指数分布族，得证。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 习题6.2\n",
    "&emsp;&emsp;写出Logistic回归模型学习的梯度下降算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**  \n",
    "对于Logistic模型：$$P(Y=1 | x)=\\frac{\\exp (w \\cdot x+b)}{1+\\exp (w \\cdot x+b)} \\\\ P(Y=0 | x)=\\frac{1}{1+\\exp (w \\cdot x+b)}\n",
    "$$对数似然函数为：$\\displaystyle L(w)=\\sum_{i=1}^N \\left[y_i (w \\cdot x_i)-\\log \\left(1+\\exp (w \\cdot x_i)\\right)\\right]$  \n",
    "似然函数求偏导，可得$\\displaystyle \\frac{\\partial L(w)}{\\partial w^{(j)}}=\\sum_{i=1}^N\\left[x_i^{(j)} \\cdot y_i-\\frac{\\exp (w \\cdot x_i) \\cdot x_i^{(j)}}{1+\\exp (w \\cdot x_i)}\\right]$  \n",
    "梯度函数为：$\\displaystyle \\nabla L(w)=\\left[\\frac{\\partial L(w)}{\\partial w^{(0)}}, \\cdots, \\frac{\\partial L(w)}{\\partial w^{(m)}}\\right]$  \n",
    "Logistic回归模型学习的梯度下降算法：  \n",
    "(1) 取初始值$x^{(0)} \\in R$，置$k=0$  \n",
    "(2) 计算$f(x^{(k)})$  \n",
    "(3) 计算梯度$g_k=g(x^{(k)})$，当$\\|g_k\\| < \\varepsilon$时，停止迭代，令$x^* = x^{(k)}$；否则，求$\\lambda_k$，使得$\\displaystyle f(x^{(k)}+\\lambda_k g_k) = \\max_{\\lambda \\geqslant 0}f(x^{(k)}+\\lambda g_k)$  \n",
    "(4) 置$x^{(k+1)}=x^{(k)}+\\lambda_k g_k$，计算$f(x^{(k+1)})$，当$\\|f(x^{(k+1)}) - f(x^{(k)})\\| < \\varepsilon$或 $\\|x^{(k+1)} - x^{(k)}\\| < \\varepsilon$时，停止迭代，令$x^* = x^{(k+1)}$  \n",
    "(5) 否则，置$k=k+1$，转(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from pylab import mpl\n",
    "\n",
    "# 图像显示中文\n",
    "mpl.rcParams['font.sans-serif'] = ['Microsoft YaHei']\n",
    "\n",
    "class LogisticRegression:\n",
    "    def __init__(self, learn_rate=0.1, max_iter=10000, tol=1e-2):\n",
    "        self.learn_rate = learn_rate  # 学习率\n",
    "        self.max_iter = max_iter  # 迭代次数\n",
    "        self.tol = tol  # 迭代停止阈值\n",
    "        self.w = None  # 权重\n",
    "\n",
    "    def preprocessing(self, X):\n",
    "        \"\"\"将原始X末尾加上一列，该列数值全部为1\"\"\"\n",
    "        row = X.shape[0]\n",
    "        y = np.ones(row).reshape(row, 1)\n",
    "        X_prepro = np.hstack((X, y))\n",
    "        return X_prepro\n",
    "\n",
    "    def sigmod(self, x):\n",
    "        return 1 / (1 + np.exp(-x))\n",
    "\n",
    "    def fit(self, X_train, y_train):\n",
    "        X = self.preprocessing(X_train)\n",
    "        y = y_train.T\n",
    "        # 初始化权重w\n",
    "        self.w = np.array([[0] * X.shape[1]], dtype=np.float)\n",
    "        k = 0\n",
    "        for loop in range(self.max_iter):\n",
    "            # 计算梯度\n",
    "            z = np.dot(X, self.w.T)\n",
    "            grad = X * (y - self.sigmod(z))\n",
    "            grad = grad.sum(axis=0)\n",
    "            # 利用梯度的绝对值作为迭代中止的条件\n",
    "            if (np.abs(grad) <= self.tol).all():\n",
    "                break\n",
    "            else:\n",
    "                # 更新权重w 梯度上升——求极大值\n",
    "                self.w += self.learn_rate * grad\n",
    "                k += 1\n",
    "        print(\"迭代次数：{}次\".format(k))\n",
    "        print(\"最终梯度：{}\".format(grad))\n",
    "        print(\"最终权重：{}\".format(self.w[0]))\n",
    "\n",
    "    def predict(self, x):\n",
    "        p = self.sigmod(np.dot(self.preprocessing(x), self.w.T))\n",
    "        print(\"Y=1的概率被估计为：{:.2%}\".format(p[0][0]))  # 调用score时，注释掉\n",
    "        p[np.where(p > 0.5)] = 1\n",
    "        p[np.where(p < 0.5)] = 0\n",
    "        return p\n",
    "\n",
    "    def score(self, X, y):\n",
    "        y_c = self.predict(X)\n",
    "        error_rate = np.sum(np.abs(y_c - y.T)) / y_c.shape[0]\n",
    "        return 1 - error_rate\n",
    "\n",
    "    def draw(self, X, y):\n",
    "        # 分离正负实例点\n",
    "        y = y[0]\n",
    "        X_po = X[np.where(y == 1)]\n",
    "        X_ne = X[np.where(y == 0)]\n",
    "        # 绘制数据集散点图\n",
    "        ax = plt.axes(projection='3d')\n",
    "        x_1 = X_po[0, :]\n",
    "        y_1 = X_po[1, :]\n",
    "        z_1 = X_po[2, :]\n",
    "        x_2 = X_ne[0, :]\n",
    "        y_2 = X_ne[1, :]\n",
    "        z_2 = X_ne[2, :]\n",
    "        ax.scatter(x_1, y_1, z_1, c=\"r\", label=\"正实例\")\n",
    "        ax.scatter(x_2, y_2, z_2, c=\"b\", label=\"负实例\")\n",
    "        ax.legend(loc='best')\n",
    "        # 绘制p=0.5的区分平面\n",
    "        x = np.linspace(-3, 3, 3)\n",
    "        y = np.linspace(-3, 3, 3)\n",
    "        x_3, y_3 = np.meshgrid(x, y)\n",
    "        a, b, c, d = self.w[0]\n",
    "        z_3 = -(a * x_3 + b * y_3 + d) / c\n",
    "        ax.plot_surface(x_3, y_3, z_3, alpha=0.5)  # 调节透明度\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "迭代次数：3232次\n",
      "最终梯度：[ 0.00144779  0.00046133  0.00490279 -0.00999848]\n",
      "最终权重：[  2.96908597   1.60115396   5.04477438 -13.43744079]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAAD0CAYAAABO8xCHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABnSUlEQVR4nO39eXhb9Zk2jt9Hu2RL8iqvsR3HWezEjp3YYQklHfY1CUuBMt1pv7QzpV+ulim8A1dfvr+2DMMLM2VogfZtOjAtFEjCUgplKy0QlpCAHS/xvu/arV1HOuf8/pA/J5KsXUfeovu6ciXRcvTR0bnP83ye5X4ojuOQRRZZrH2IVnoBWWSRhTDIkjmLLNYJsmTOIot1giyZs8hinSBL5iyyWCfIkjmLLNYJJHGez+atssgi86CEOEjWMmeRxTpBlsxZZLFOkCVzFlmsE2TJnEUW6wTxAmBZZCEIfD4fpqam4PF4VnopKwaFQoHKykpIpdKMHJ+K02iRjWZnIQhGR0ehVqtRWFgIihIkeLumwHEcTCYT7HY7Nm7cGP50NpqdxdqBx+M5a4kMABRFobCwMKOeSZbMWSwbzlYiE2T6+2fJnEUWUTA/Pw+z2Rz1eZZll3E18ZENgGVxVuCFF17AQw89FPM1d9xxB77+9a/z///73/+Od955B//3//7fJa89efIk/ud//geff/45enp6sH37drS0tOCxxx7Dc889B4/Hg2984xtCf42YyJI5i7MCN910E2666aaQx7Zt24a+vj7+/7/85S+xfft2qNVq/jGfz4dzzz0XAOB0OvGP//iPuOeee3D06FF85StfwZ133omf//znOHToEJxOJ/7P//k/2LBhw/J8qTBkyZzF6oXBAIyNATU1QHFx2oezWq2wWCyRoskAALlcjocffhizs7PYvHkzvvCFL+Ctt95CWVkZGhsb8d5776G9vR1+vx/PPfcc7r33Xhw5cgSDg4O46667cODAAfT09KwYmbN75ixWJ/74R6C6Grj00sDff/xj2of85JNP8Oijj0Z9/jvf+Q4uueQSXH311Xj88cfx2GOP4dixYxgdHQUA7Nu3D3feeSeefvppmEwmAMDLL7+MX/ziF6BpOub+ejmQJXMWqw8GA3DbbYDbDSwsBP6+7bbA4wJidnYWra2t/J/Z2VkcPnwY3/jGN3DnnXfie9/7Hux2e4jbDQT20vv378fQ0BD8fj9yc3PR39+PlpYWQdeXLLJudharD2NjgEwWIDGBVBp4XAB3m6CsrAwnT54MeezWW29Fc3MzPvnkE5xzzjkYHR1FVVVVyGsef/xx3HHHHXjrrbdwzz334Ic//CFUKhU0Go1ga0sFWTKvADiOg8/ng0gkglgsPuvzr0tQUwPQdOhjPl/g8Qzi/vvvxxtvvMH//+GHH8bMzAz+8R//kX9s165dePzxxwEA//RP/4Tc3FzU1NTgwIEDGV1bIsiSeZnBsixomg6pBBKLxZBKpZBIJFlyAwHre+hQwLWWSgNEPnRIUKscCffffz/uv/9+AMDp06dx2223Qa1W45prrsF3v/tdFBUVLXnPoUOH4PF4cOmllwIAnnrqKfzhD3+ASLT8O9gsmZcJHMfB7/fD7/eDoiiIxWL+cZZl4Xa7eRJ7vV6o1WrIZLKzl9xf/jJwySWCRbP379+PmZkZAEBraysAIDc3l/83AGzYsAH79+/Hiy++CK/Xi8cffxxNTU144YUXcNVVV6G5uRk/+tGPsHXrVgCBevM33ngDzzzzDADAbDbjwIED8Hq9OHToUFrrTQXZRotlAMdxoGkaLMvyxPT5fFFf293djZqaGqhUKgDrw3L39vaivr5+pZcRE8888wwmJydx8OBBbNu2LeQ5juPwyiuvQKFQ4Iorrkj5M6KcB0F+0KxlzjD8fj+mpqbAMAwqKipAURQ4jgPHcRFJSVEUb7nFYnFEyy2RSPg/a5XcqxHBe+NwUBSFgwcPLt9iUkCWzBlCsFvNsizvXicLQm6yB+M4DgzDwO/386+RSCS85RaJRFlyn6XIkjkDYFkWPp+Pd6uJNU4UsV5PjkcQTm6KokIsd5bcZw+yZBYQhFhkP0ysaTRyRnO1k0Ekcvv9fpw8eRI7d+7MkvssQpbMAoHkjhmGWUKwcDLHs9bJWvLw91IUBYZh+D233+/nbzBZcq9fZMksAEjumFjacHKkQ850Ecly+3y+JeSWSqUQi8VZcq9hZMmcBsJzx9EKBYTcM6eL4Bw3EJncJJgmkUgi3pzWOhLd3szPz0MqlaKgoCDi8yzLrkhxSDSsnpWsMZDcMSFyrIsjFXIulyUPToORNBdN03A6nejt7YXJZILL5eIDeivlYQiF999/H1deeWVCr/373/+Ou+++O+JzJ0+exJ133okLLrgA+fn5uOCCC3DHHXcAAJ577jk89dRTQi05YWQtcwogQa5obnU4UrHMK4Vgy22xWKDT6UDTNLxeL+99hLvlawUnT57EN7/5TZSVleGCCy4AAHR0dOCPf/wjrr32Wvzyl7/EE088kRUnOBuQqFsdjtXkZicL0gwCnPEWaJoGvdgIkUlyC6lN8Pzzz+Pee+/FK6+8gh07dgAAfvWrX6GtrQ3XXnstgLUvTpAlc4KIlDtOFKuJnMkgfM3kOy8Huf/4x0CfhUwWaKA6dChQrp0qcnNzcemll+KWW25BXl5eyHN79+6FVCrF3//+d/h8PpjNZtx5553o6OiAwWBAa2srGhsbsW/fPuzbtw+HDh1aIk7w1FNPZcUJVjuINfZ6vSkRGVjbljleLCB4zy0SifhYgtPphM1mg8lkgtVq5UtYE0EmtAmuvvpqlJSU4OGHH8axY8fwyCOP4Fvf+haOHTuGv/71r/zr1rI4QZbMMRAc6SVudaolmauFnMkg2TWTcxRMbrPZjImJCb6ghmEYsCwbM5hGtAmCQbQJhILFYsHU1NSSx2+99VY88sgj6OnpgUQiiSpOIJFIsuIEawUsy8Jut6O/vx9NTU1pBaXWq2VO9L3ELQ9uMgl+DXkd+fdyaBOQ5pdgZMUJ1hmCSzJJmkaIksvVQs5kIMSaScSfHCv8XEYid2EhhUOHqEVtAkowbYKDBw+iuLgY09PT+NWvfoXa2lpcc801uPfee/Gzn/0MF1xwQVacYL0gvO+Y7AHTRSQyZ6qcU2ikeyMLJ3O84xNyf+lLwBe/CIyPU9i4kYJOR4HjUl/P3NwcBgcH8cgjj6Cvrw+PPvooLrzwQnR2duJf/uVfoNVqUVZWhvfffz8rTrDWEakkk2EYnDhxgs8zpgqHw4Hh4WHs3LmTf4zIBkW6gw8NDaGgoCBq5VEiOHHiBNra2lJ+PwB89tlnaG5uDqkYSxbT09MAArn58Ib/RBAcMZdKpXzcItg1TwQvvfQS2tvbccUVV+D8889f8vzRo0fR3d0NuVyeFSdYq4iVOxbKQoYfZ3Z2FoODg6AoCrm5ucjPz0dBQQEUCoWgnysEhLDM6bic5PMZhoFMJou6545H7uuuuw7XXXdd1M+54YYbcMMNN8RcR1acYBUjXu5YaDIzDIO+vj7QNI22tjZQFAWn0wmLxcI/rtFoQNM0cnNz0/7cdCHknlkoJLLnTsVyrweclWQO7zuOljsW6kKgKAo+nw8nTpxAWVkZqqqqeAUStVoNtVqNqqoqsCwLm82G4eFhjI2NYXJyEnl5ecjPz0deXl5a7m46a08HQpE52T138PNnC7nPOjKHu9XL8eMajUbYbDa0tbVBq9VGfZ1IJOLJq9FokJeXB6vVCrPZjJGREYjFYuTn5/PPZzpiKoRlDu4sSpfYibw3Hrk5jgvRTVtOcmd663RWkTmdksxUP6+/vx9OpxMajSYmkYNB3HKJRIKioiI+JULTNCwWC+bm5jAwMAC5XM6TOzc3NyPfRyjLLJPJYDabUVBQsKwECv8sl8vFq56S55fDcnMcB5PJxMdFMoGzgszhuePlyAG6XC50dnaipKQEtbW16OzsTPi90fbqMpkMJSUlKCkpAQC43W5YLBZMTEzA4XAgJyeHJ/dqCaARMhcVFcFoNMKQYk0mTdOQhZeFZeg4mSK2QqFAZWWloMcMxroncyw5n0xBr9djcHAQDQ0NyM/P59slhYZSqYRSqUR5eTk4juODaUNDQ3C5XOjt7eUj5UIQIRUQMkskEpSWlqZ8HCFSbQDw6aefYs+ePRGfI7LGBGTda0WFZV2TOZ6cTyY+b2BgAE6nE21tbTyBlqOck6S5cnNzsWHDBnz66acoLy+HxWJBd3c3GIaBVqtFQUEB8vLyIJEsz08vdDQ7kwju5QZiq7CsRnKvSzKn2necDtxuNzo7O1FcXIytW7eG/MgrkTemKAparRZarRY1NTVgGAYLCwswm80YGxsDRVG8S67VajN2jtYSmcMRidxEqGFmZgYlJSVQqVSrRmJp3ZGZ4zjMzc0hJycHMplsWU6uwWDAwMAA6uvrI1ZtRSKzSCQKcenivT5diMXikKoyn88Hq9XKbwmI1lV+fj7UarVg5201kTl4PFAqCCa32WxGSUlJiAoLsdwrpcKyrshMJGWnp6exYcMGyOVyQY4b7YJkWRZDQ0Ow2WxobW2N+nmRyLnSASqpVIri4mIUL3YveL1emM1mTE1NwW63Q6VS8X3JKpUqZRKsNjILRTAiZRxNqOH73/8+7rvvvpRKWFPFuiBzuFtN9KKFALGg4QUbHo8HnZ2dKCwsxO7duwWvclpussvlcpSVlaGsrAwcx8HlcqGjowMjIyNwuVxQq9V8MC2Zm6QQZBbqXAhJ5vBrIrzd02w2Q6lUCvJZiWLNkzlS7jiWC5ssIhHLaDSiv78f27ZtQ2FhYULHWEugKIrfpjQ2NoLjONjtdlgsFpw+fRo+nw9arZbfc0ul0qjHEoLMQpEw0k05HcT6Xk6nc9lLctcsmcNLMoN/bLFYLBiZg28MHMdhaGgIVqs1pludKKJd6Kup0QIIrEej0UCj0aC6uhosy2JhYQEWiwWTk5PgOA55eXkoKCiAVqtdEjRaLSWhQlrmeL+Py+XKkjkRxMsdi0SiJSoSqYKQ2ev1orOzE3l5eWhtbV1z1lZIiEQi3ioDgViF1WqF0WjE8PAwX3ZaUFCQdtAJENYyL1dQyufzLXtuf82ROZHcsdButsViwcjICLZu3RpRbUJorDbLHA+Ryk7NZjNmZmZgNpvh8/lQVFSE/Px85OTkJE3udNsoCYQk82q8ma8ZMieTOxaKzCQQNDY2ht27d2e0rnY9QSaTobS0FKWlpfD7/aisrOTPI9lLEsueSJBICOsOnIlAp4t4N9qViuCvCTKHy/nEO1FC7JlpmkZnZyc4jsOOHTtSJrJpwQG5NLnTTFGUYJ7FSoPjOKhUKuTn56OiogIcx8HhcMBisWBgYABerxcajYYndyTXdLW52YncFFaC0KuezMmOggEClpkExlIBidpu2bIFc3NzSbu8Xp8f/eOz6BqewqzRCrlMgjyRF+cl+AOvRhcuVYRf1BRFRezhtlgsmJ6eBsMwIT3cEolk1QXA4pHZ7/evSO/5qiVzOiWZqbrZHMdhdHQUBoMBu3btglKphF6vT5jMk/MmdA9PYWByHj7/mQCcl/bj5OQMZH/9FJfu2YECTfx941raM8dCPBeZ9HDn5eVh48aNYBgGVqsVFouFLztVqVR8+jEdMgpJ5ljHcTqdyMnJSftzksWqJHO6fcepkJmmaXR1dSEnJwdtbW38jxXP5bW7POgZmUL3yDSsdlfMz5icN+N/Xj+Gc3dswrYNxfD7fBFLJ9daACwWkrWqYrEYhYWFfP6eVPTZbDacPHkSMpmMLztNtoc7HgkTRbx89UrkmIFVRuZE5XziIdk9s9VqRU9PD+rq6vheYYJINwaGYTE0NY/ukSmMz5nAsokTz8+weOPDz/E67cYFjbVQipFyddVaQLouslQqhVarhc/nw+bNm+HxePgpGQ6HAyqViie3UqmM+VlCFY3Ec7PPesssZN9xonlmjuMwPj6Oubk5tLS0hChQEARbSYPFhq7hKfSOzcLtpZe8NpHPM5mM8Pv9KC4uRvuEFS1bqrG5rBS2BStOnz4Nv98PmUwGmUwmWPR1JSF0BZhCoUB5eTnfw+1yuWA2mzE0NASPx8PfGPPz85fcGFmWFaT1M97v4nA4zl7LzLIsDAYDNBqNIG1kibjZPp+P10nes2dPVPfLx7DoHJ7CzMkBzJkWUl4Tx3GYmZlBTo4KRUVFi14Ii/aBCQxPG3Bxaz1aWgKtiqOjo1hYWMDnn38OiUTCdztlShook8hkBRgpO83JycGGDRvAsiwcDgfMZjN/YwwuO12uPTNRfVlurCiZg93qnp4e7N27V5DjxiPzwsICuru7sWnTpojqFxzHYWIxmHX8VC/EEglyclK/07rdbvh8PhQXF0GpJNb/jGtuc7rx0nufo76mDPtatkKj0UAikaCmpobvZgqWBiLkXit57+WqABOJRHzZaXAPN5FW8ng80Gq1kEgkS8pOk0E8y7wSpZzACpI52dxxMoi2Z+Y4DhMTE5iZmUFzc/OSu+eCw4WekWn0jE5jweEGADAsB1HKwSgOVusCnE4HpFJp3AKJ3rFZjM0asaO6GJUFgYshvJvJ6XTCbDajr6+Pb3hYbvWQ5UaqFjW8h3tgYABSqZQvO5VIJHysIpke7kQCYGeNZc60nE+kPbPf70d3dzekUin27NnD/xh+hsHg5Dy6h6cwMW9COG9TjSyzLAu9Xg+xWIzy8oqI40Mjwe314b2OQRSplbiluATanDM3gGBpoKqqqiXqIaRmurCwcN1EwwFhtbdJMwgQ6OEmo10dDgcUCgUfTIvVw80wTMwb51kRzV4uOZ9wN9tms6G7uxs1NTUoLy8HAMyZFtA9MoXesRl4aX/UY1EUlhA8Hmiaxvz8PPLytFCrNSHHSfSanDEt4A9/+QjnNm5Cy+aqiOcq3PKQmujp6Wm4XC50dXXxzy93b62QyFQFmFwu58tOOY7j1U5JD3dubi5P7uAtDcMwMbMOTqczrTlhqWLZyJxI7lioOzAhM8dxmJqawtTUFJqamiCSyPBZ3xh6Rqfg8fpgc7oRb2ZXspbZ6QwEYHS6krAfnEJgn5zI9wtMPKT9DN5vH0D/+BwuaWuALj/2MO/gmmi73Y7a2lqYzWa+bJK45Pn5+WvKJRcqCh3rpkCKU1QqVUjZKdnS0DTNB9PijQ9yOp3YsGFD2utNFhn/RRPNHUdT9EgFxM3u7OwERYlQXLkR77YPYnhKDybIYqsUcrAsC08MyxwgX/ycNcdxMJtNoGkaFRUVEIlCv0eyljkY82Yb/vj2cezeWo1zd2yCJIFzFCnSS1zy8fFxXtCvoKAg4nSM1eSmr0Q5Z3DZaXgPt8FggNlshtVqjTg6aF262cmMghGLxYLlVZ1OJ+aMZlAuQG/3wu6KvF91eQKD1DU5CtidXnARJtgGKsBiX9gM48f8/DyUSiXKysoQ2foSyxz8/8gInKbQz2RZDid6xzA0pcfFbQ3YoEvOjQvvQfb5fCHTMch+cTW65Kuh0SL4/NE0jZKSEjAMs2R0EMuysNvtgpKZoqjXAcxwHPftWK/LGJmTLckkZE4HPj+DD0504KOOXuitTlRVxb8xcBwHm9MDhUwKSkTB7QktBonnZns8HhgMehQWFkKlih7BTG7vHf1cWewuHP3bSWzfWIEvNG+BQhZdsicWpFIpdDoddDodv18ML74gAomxZIGWA8vZ7ZTocWQyGXJycpaMDvrtb3+L1157DT09Pbjhhhvw9a9/Pa39M0VRlwNoBjAT77WCkznVUTDpkHnaYEHn4ASOfdYJ2seguLgY1ELsOulweGgfKArQ5CjhdHvALFrj6CTksLBgg91uR2lpWQIXfLhljo1YxOc4oHtkGqMzBvzD7nps3lAS/cUJIHi/WFlZybuUJpOJbwMlLnkmNbajYTV2TYUfh4wOuvfeezE4OIjbb78dExMT8PtjbeFig6KoHAD/H4CHAOyI93pByZxOSWayZHa6vegZnUbPyDRmDWbo9fPQaLQoKVAj1UH0HBco4JBJxVDIJXC6vVFkclno9QZQFFBRUQ6Kin+BJGOZI7nZkeD00Pjzh6ewqVKHf9i1DWqVMEUkpPhCqVRi9+7d8Pv9sFgsvMa2XC7nXfJ0ZHgThVAkFPKmEC/PXF9fj8suuyzdj3oUwH8AWFpnHAGCkZnjOHi93pRzx4mQmWVZDE8b0D0yhdEZA1g2oBpptVpRUqKDTCZMkwLtY0D7GKhzFKC9bgQTy+fzYX5+brHSKLGpjgEkZ5mTwfCUHpPzZlzQtBlNdZWCk0sikYRobBOXPFiGl5A7E7pXQskGAcL0iidSAaZWq9P6DIqivg6A4zjuBYqivpHIewQjMyFwqicrFplNC47FBocZON1eAAHraDAYwbLsYvRYeNfP7vRAJBJDKQucJpfLCZPJBJ1OB7k8OSuY7J456dy2z493P+tF/8QsLm5tSO7NERDLiimVSlRUVPApHLvdDrPZjO7ubrAsG6LUKQSEkg0SCvE8BYEqwO4AkEdRVB8ALQAlRVEijuO+Fe0NgrrZ6WhvhZM5XK0jGD5foChDrdZAq9UguludaF43OvwMC6eHhkhsh8PuRHl5RYpBlMxZ5mBMG6x49q1PoBZ5sWsXC7E4s/vbYBnempqaEKXOoaEhuN1uTExMpDVDejlVNRNBPHc9XoVYgp/RSv69aJkvWLFodrIgZJ7Sm9E1PIXByTnQvqWW2uGww2KxQqcrjmkdyV5XiCJ/p9MFkUiMsrIyqJRy2J2epI+TiT1zNPgZFj2Tejz71ie4pK0BZUV5SR8j1XMXrNTJsixOnjwJiUTCN4qQqqpkerdX04gbILarvpL5+VVBZrvLg86RGQxM6gFJ5B+Y4zgYjUYwjB8VFeVLijLCIRIRMqe+Lq/XC6PRAKlUisLCQjAsB7vTgxylDH4/A2+Em010LI9lDoZxwYHn//opdtZtwPlNm5MWFkwXHMct1qaf6T8Ob1EkLnl44UUwVptljgeh+w04jnsKwFPxXifor5vKF5g2WHD0bydgMlsXxdyWkjkQdJqHWp0LrbYIibjOpNgj1WvAZrPBZltAUVER7HZ7yHNONw2RiIJGpYDd7UnI4kayzNFPV/J75mjgOKBjcBLDMwZctLseteXFCb5P+D7k8KoqovdFgmnRupiECIAJaTETkdpdCay4Za4ozsdXrtiLw299iIk5w5LnSa1zcbEuqf7dVLudAh6AASzLory8AgzDRDwOy3KwuTxQyqXgAHi88dRAl98yB8Pu9OCV99uxtaoUX9y1DSpF/KhzpsfKhOt9BXcx2e12vneb1CykAyGryGKthabpFSuyWXEyA0CBJgfXXdiMTzr7MWZ2g/YxvMSOz+ePWOscD6mQ2e/3YW6OeABaANRiw0b097i9vsWSUCXsTk/EktDAes5YZvLdXC4XFAoFlEoVFApFkIggkCni90/MYXzOhC80b8GO2oqorxPCuiRr3cO7mEjvttPpREdHR4hLnmyAabmGz61ULzOwCtxsAolEgs0VRdh3bjVeO9aOE52ng8rlkj8uRYmSuiDdbheMRiOKi4uhUAT3ECfmVtmcbihkUohEgMsTyUoHLDPDMJifn4NCoUR5eTm8Xi9cLjesVgtEIhGUSuWiVnTCS08aHtqHtz/tQf/4LC5pa4A2N3JNQqYtcywE924bDAY0NTXxKbDg3m3ikscj6nKK+a1EkwWwSiwzcCaa7XU5UK0Gai7bi8+HZhJwXyMjYJkTSZNxsFiscLlcKC8vh1gcekqSsfBnSkIVcLq9fElo4DgATftgMBhQUFAIlUoFhmGgUCj5mwfD+OF2u2G3O+D1emEwGKBSKaFUKpP2TBLBxLwZv3/jY5y7vRa7tlaHEGIlLHM0kBbI8N5ti8WCmZkZ2Gw2XqUzWqOIUDK7q1WZE1hFZCYD2rxeL9ra2iCTydC0pRbvnDiNoan5lI4X74JkWQbz83pIpVKUl5dHufCSr6m2OT2QScVQyiVwLBa5+Hw+OJ0ulJWVQSaTRYy0i8US5OaqIZPJsbBghUajhtvths1mAwAoFEqoVMrFlE4iJImfZ/f5GXxwahD9E3O4tG07dAWaxe+xekaxAku9BFILXVJSEqLSGdy7HTw/er0L4AOrxM32eDzo6ekBAOzatYs/Tq5KgYP7dqFvbAbvftYLlydxedt4ZKZpL+bn9XwxQ6rHiXr8xZLQXJUcU9OzoGkfCgsL+XLHWKeKPCeXKyCXK5CXlw+WZeB2e2C322E0GiGVyqBUBsgd7k2kAr3Fjj++cxwtW6pw3o66tI8HLF9+OFrvNpkfDUDQqRhZNzsKDAYDBgYGUFdXh4mJiYg//raaclSVFuFvn51G79hsQscleeZIOFPPXRK3ljgV2SAClmUwODwKlVIJXVE+uIQv7EjiDWL+ggU40LQPbrcbBkMg8q5QBNxxhUKRMoFYlsNnfeMYnjbg/O01q8oyJ4NIvdtTU1P8VIxEtb4iIRHN7HVhmZMBy7IYGhqCzWZDa2srJBIJRkdHo75epZDh6r3N2Fpdhnc+7eHd12iIJCoQLEKfeD13ahcj0QEjlt9sNqEgTwmpRASfPxHlkthrIkL5Wq0WHMfC7fbwrqZEIubbUANpkuS+g9XuwisfnEK+gkLd5q1QyFNLtayWyi2pVAqNRgOWZVFbWxvSKOJ2u/lGkWhTKIORDYCFwePxoLOzEwUFBdi9ezfvyiZS111XWYJKXQHea+9D11B0xctw99jv92N+fi6tCHmiICmV0E4uCk63FxqpHJocxaL+WLS1A8ns0ylKxPcjA4EU2+zsLMxmC/x+PxQK+aLVVibhYnIYnbPg6b98iC+2bMPW6qX64vGwmiq3yFoi9W6TKPnU1FTc3u2zhsyJ3IVNJhP6+vqwbds2vlgg0fcSKGRSXH5OI7ZWleHtT7t5jevwtRAyu91uGI0GFBUFi9BnAhwsFgvcbs9iZPzMj07cdYZhYXN6oFLI4aXpiPXn6UIikUIsFqOkRAeOA7xez2L6awEiEQWlMuCSB6xQ7PPu8tB4/eNO9I3P4h92b4MmJ3FJodVimYHoNxaRSAStVgutVouNGzdG7N0m5M7JyUkoAKbT6TL5VaJi2Swzx3EYHh6GxWJBa2urIAPSasqK8PWrLsAHHf3oGJwIcU0pSgSWDXTwOJ0OlJWVZ1SRMqCTPQ+JRIqysrIIF3FoVNzl8YJjOahVcjjc3jC3WrhyToqiIqa/FhYWQNM0ZDJ5xPRX4PPPfIexOSNeeu9z7Ny8ATvrNiRE0rVA5nBE6t0m42WdTidEIhHUavXiuVvqkq+baHY0eL1edHZ2Ii8vD62trYL+wDKpBBe3bcfW6jK8ebwbFptz8RkONpsdSqUS5eUVGb2oiGCBVpsXtSk9UiCN5VjYXd5ASSjHxVEJFQYk/ZWbqwbAweul4Xa7YLPZwHHgI+TkfEnEIqgUMrg8NMw2J/72WR/6xmZxSdt2FOXFdieFSm8JgVRdfqVEAmVxMd8o0t/fzw9UYBiGD7SRRpF17WabzWb09vZi69atvPhZJlCpK8DXr9qLDzsH8XHnAMxmC2QyGX+HzRRcLhdMJmMCggXR89WBktBAsYnD5cloOWf4muRyOeRy+ZL0l99HQyIWYcFmg8erCPFqZk0LePatj9FavxF7GjZGlf5dTbnqpLW3OQ6izz6D+NQpAACzbRvYc8+FVCpFcXExCgsL+d5tk8mEkZERvP322+jt7cX27dvTihdQFCUD8AsAlyBw4dzDcdzReO/LmGXmOA4jIyMwmUzYvXt3wk0S6ZwEiViMrWV5sOtzIBeXQ29eSOk4iYFbdOEjV46FI16KixSbyGUSSDIsKBANIpEYhflaSIoKYF6wwWKxwu9nYLeT9BepI5eDYUU43jOCwcl5XNLWgIri/CXHW01kjjeFIhzUyAjEn30GtrISACDu6QE0GjByOX99BvduA4BOp0N7ezsOHz6MRx99FEePHsXmzZtTWW4BgHc5jvsniqK2APiUoqg/cRwXsxxScDJTFAWv14uuri7k5uaitbU1KYXOVMnMsiwGBwfhdDpx+UUX4lybDW9/3I45OxMifC8EAoJ+eohE4hiVY+FIrJLMS/vBMgyUcinEFAVmmdrpcha7qJweGoAfHBf4PUhwKFL6S6lUwufz4fC7J9C4qRIX7NwS0jMt9Gzm5TwONT8PLjcXpIeW1WqB2VkwGzZEjWZXV1dDpVLhgQcewJYtW1JeK8dxcwCOLP57gKIoPwKifjGtk+Bktlgs6OnpwebNm5OO6pH67GQDVWRPXlBQgJaWFlAUBalEgqaNZbikfAPePN6d1mzlAAKlkaSzikjlJIpkik84BJo1isQiKCUSODyxc+qpgqKAHIUcDMsuknjp82f+vTT9RfK1fr8fRqMRXQNjuHJvM7ZWB+Z5CWWZV0QAX6MBXC5gsRaccjjA1dbGrQBzOBwJNX4kCoqivgmgk+O4uBew4FK7c3Nz2LVrV0pTEVLRzrZarejp6cGWLVtC9sdkRE1xvga3XnYuTvaN4aPOQfiZ5K00IaLHE0hxhXdWJXgUJDLmhnwewIH2M6D9DNQqOdxeX0prjwSRiEKuQg6vzxej+Cb2nUcikUKtlkKt1iwqs3pgXrDht0ffQpUuDxe3NSBHkX7GQigxv2TJzG7dCmpyEqLpaQAAV1wMtrERzMBAxpU5CSiKugfAzQCuSuT1ggfAGhoaBBP1iwWO4zA5OYmZmRm0tLTwFoMgWFxQJBJhT0Mt6ipL8OYnXZg2WJJcGYWFBSucTmfKKS6KAlL19u0uLyRiMdRKOexxKt9iITgybXMlomOWGImC01/5+YCX8ePVT3pRW6hCWZ4Sbrc7ZSneFZtmIZWCuewysGYzwHHgCgoAiSQhmd3wazEVUBT1KwA5APZyHJfQRIcVr80ORqJkZhgGPT09EIlEaGtri3hyIymFFmhycMul56BjYAIfnOpPqGAjMC/LB6/Xm2aKK5nuq6V5Zj/DwO5mkKuUg/b5QCdQEkogk0ggl0ng9HhhS1CMMJ2tulgsgVgsQd+cHSY3g6tq6+B2u9HV1QUyIznawLql61jBaRZiMbiwbEi8ohEh+qYpijoXwFaO4y5J5n1rjswulwunTp1CZWVlzLGZJJgWDoqi0LK1GrUVxXjreDfG50xRj0FKQEUiEQoLi9K6qJJ5a6zXOtxeiEWioJLQ6C9WyKSQSsRwuL2gXcnmsNMTQyTHmLfY8cqH3Thney1aW1rAMEzIwDqlUonCwkIUFBREzHishqFxiR5HQO2vZgCtFEUNBT32fY7j3oj1poxEs1NFPDKTErsdO3bEFVgne+Zo0Oaq8KWL96BzaBLvtfctGbhOBsIVFRXDarUi/bxvsq2U0V/LsIGSUKVCDpZh4fWFrj1HIYNCJoWH9i1LIUo0kCoyP8Piw84hDEzO49K2BpQEDawL5OkDJb4+nw95eXkoLCyEVquFWCxeuQBYFMRT3hRCmZPjuCcBPJns+9aEZeY4DkNDQ1hYWOCFC+IhUUH+proNqC0vxtsnejA8pQdwRpmzrKwMEokUCwvWtMsrw39fm80Oi8UcVE6pCnLPErsY3B4aItGZYhNVUGTaQ6em0EIQXs6Z2jFCXWSDxY7n3vkUzZurcN6OTZBJJXxbZ1VVFa/WSQT05fJAg4jf70/b3V6Opg+O47K62QSRyEzTNLq6uqDRaPgOq0SQjKhArkqB6/btRs/IFF544324PV5UVFSADIRLVaAgbEWLx+BgMpng8/lQXl4Bv98Hl8sNvV4PgOO7mxL/PArgAKVCDj/LpiyzlBksJSDLcvi8fzwwZ7q1HjVlZ6oCw9U63W43JicnYbPZcOLECWi1Wr5VcaUE/WILXtCC9BykilXtZi8sLKC7uzulnHWy6/B6vXAaZ3DzRa0YNjrRPz4XfDSk62YTYcDZ2TnI5TKUlJSCZRnIZHLIZHLk5eUtllO6Ybfbeb2waBpgYpEIOUoZXB4fH5mmKECtUsDpTn7ixlKkv2eOZd1tTjdeeu9z1NeUYV/LVijlS70tsp+WSCSoqanBwsICzGYzxsfHIRKJUFBQgMLCwoTG3iyHZV7JJgtgFVpmny9gWaampjA5OYnm5uaMnyBy0yD14w0AtlXP468nAiIIQlhmv5+BzWZHcXER3+QQjoCaSC5ycnJA0zQ0Gg3fBAEASqUKWk0uNLk5cHroJZFpjgtMB5HLJFAs8/SKyIh/Q+gdm8XYrBH7WraivqZ8yfMkzxyuHkLTNMxmMz/2JpFJlJmuRltJlRFglZKZTBPcs2ePIPKosTAzM4Px8fEluerNG0qwoaQAf/+8D38zGNIis8fjhsVihkqlWiRyPFCgKIQ0QUjFFGivF0aTBTOzc4v7SdWi1Q69wLy0Hx6fP2ClPfQSxZVEIMyeObFjuL0+vPFJN/rGZ3FRawO0QT3T0Qgkk8lCNLbDJ1Emk/5KFKtZmABYZW62z+fD9PQ0amtrUVVVldG2RZZlMTAwALfbjba2toh7MIVMiivObYTU50TXpAn+FPgcCKbZUFhYCI8n1JKSAF3w3jwcKoUMFBVQKaEkpAuM4/W2FxYWFtUzAkE0meyMTJDd5YFMIoEsSCV0eZGcqz42a8If/vIRzmusQ8uWKt4jikfGSJMozWYz5ubmMP7RR1C73WC0Wng8nqSmooRjNQvgA6vIMhuNRgwPDyMvLw/V1dUZ/SyaptHZ2Yn8/Hxs3bo17k2jUpePhs0b0TNpxKkwEYTo4GA0mhb1xgJi98HPEYsTiIAGpmZQ1BlCK2RSKOXSKIqkFK/cmZ+fzwsOWK1W+Hw+yOVysCwLlmXTKAnN7J45Gmg/g/fa+9E/MYdL2hpSKueUSCTQ6XSo/M1vIHv0UXBSKRgA3f/xH7Bs2pTQsLpISERlZF1Z5mQR3CrZ0NCA2dnE1DeTOX7wxWC329HV1YW6urqEg2oikQgSsQiXtG3HtiUiCEvBsizm5+cglytQWloCckGfIW6gaUMkCn4coMAFyi3dXri9dMAaJ5C3DBYcCNRJe+FyuTA3N7c4JUMFn08JpUKBXKV8Ga106jeEOdMCnn3rE1QX5mDXlqqk3y/+5BPIfvlLUF4vKK8XIgA7f/IT2Hp7Q4bVkQmfBQUFcZU6V7MyJ7DCbrbP50NXVxdUKhVaW1vhdruTbrSIBZJrJj/A3NwcRkZG0NTUlNQdNKD0GbBolboCfO3KgAjC5/1jS/ajRHUkL2+pHveZPCQVcpFLxGLkKOVweWgsOFzQ6/XIzy9YfM+ZWVeBcxu/aEGhUEAiCbRn+v0M3G5XQNfK54dCoQj0LEul8DPRXYzl3DNHA8tyODU8jQm9FdddrMIGXUHC7xX19S15jJqfh5hlQ9NfdjvMViuv1KnRaFBYWBgx/XXW7ZkTBbGQtbW1KC0NKD+m0jUVC4TMIpEIg4ODsNvtaGtrS3pKH3GHCaQSMb64axu2VpXizeNdMFodAM4IBy5VHeEgkUhA0zRmZ2f5VsIclRJKuRQOFw2b0w2Xyw2z2QSdTsdHZInVJlYd4HiSJGK1JRIJ1GoN393k8Xhgti6A9nqRq5SDFUmgUqkgkWRicqEQrnpg2ubRv53EjtpKXLBzMxSy+GtlI/QTczodQH57mobkhReQd/w48iQSVB08CN8XvgCbzQaTyRQx/ZVIk8VZR+aZmRmMjY0tsZBCk1ksFvP7Y7VaHTItIxkEW+ZglBXl4atX7MUnPcN4++N2WK0LEbqquMU50SJ+RCzrp+FxOqDX6yGXK/hpC4GurLIwVU9ixSkAokVJ4kDxCcexfCdWwGWPb7WJMicQ6EnmGD/sC1a4PDSUyjMTKVdqz7z0GNxiIAzoGp7C6IwBX9y1DZs3lMR8H3P++aC/8x3Ifv1rcIvjabzPPss/L37jDYg//hhcVRXg90Py3HPgdDrk1dcjLy8PwNL0l1QqhVQqjSrm53A4UF6+NL22XFhWMrMsi/7+fng8HuzZs2eJGyM0mVmWRXt7O+rq6njrnwpilYZSFFAg9WPf9g2YWCiDYdFKB8AFRaypoMi0GDKFEuq8fLjdHpjNJvh8AX1rh8MJlUoFaZQ8MUVREItD99pn9uFnrHYiQTqJRApIpChUKrFBLofBbOGVRCgqYNX9fn8aqqbC3xAcbi/+/OEp1FXq8A+7tyFXGT06Tf/0p/B9+9vwT0+jF8COtjb+OfHp0+CKigJKIjIZoFCAGhkB6uv514Snv8bGxviaBJZlkZ+fj8LCQj795XK5zo49s8fjwalTp6DT6bBt27aIrwt3Z9OBXq/HwsICGhsbUVIS+y4eD9GKRnw+Hzo7O6HVavHFveeC4zic6B3FR11DfD0xAKhVSjAsuyQyzXEcFhYWoFLlID8/f7G00wWj0QCGYReVMgOaW5HO1xmrLeaPx7IcaJoGx7FgGDYoQh7darNsgCT5eVrkazXw0H4sLFjh8XhhNBrBskzQ0LrEx98EIvSZ0QAbmtJjUm/GBU2b0bipMurncNXV8BYXAyMjoY8XFUF0+nRAGggAvF4gf6mOGQFFUbxAZEVFBa+vTbq/Tp48idOnT6OuLv05XRRF3QTg3wEwAB7gOO53ibwvI5Y5/OInCp319fX8SM5MgUTHzWbzouh98oon4YjUgeVyudDR0RGy56coCnsaalFbXoS3PumCzeWFz8/AEaG8MtBeOQ+NRgu1OnBBSaVSXnOLZTm43S44HIFBcTKZlN9rR9u3URQFmg4QsKSkZLHrKNRqxyK2J0gl1OmUQKmkoNGc0f9yOp3Q6y2QycTIyVHFHVpHgn3pIbp199J+/PVkL/rGA2msAk1kqxip8MR/4ACkY2OgJiYAjgO7bRuY1taYK2EYho+3BOtrcxwHtVqNd999F08++SQee+wxPPfccymJ+VEUpQbwCIBzESBzB0VRr3IcZ4j33oy62RzHYXx8HPPz80kpdKYKv9+Prq4uKJVK7N69G319fSmrngRDJBLB7z/TSkha9sJbMcmInbxcFW665Bx0DE7gw86hJcfzeAKzl4uLi6KeE5HozGRDjgvs3wLppsB4W5UqYLVlMjl/sTscDiwsLKCsrDTINQ612gGCBUfIQwtWiEqoVCyGVCziXyOTqTAykguTiQLHMaiocEKtNoBlOX6vLZeHexDCBMDiWfdpgwXPvPkx2ho2om3bRojD1E0jkZkrLgZ9zz0QTU6Ck0jAbdwIxNlOBOZpL/29KIpCfX09ioqK8G//9m+or69POsgahMsBvMdx3PTisd8FcDGA5+K9MWNkJkLhMpkMbW1tGS9yJ6IF1dXVfBAi0TbIeAgOgBGpovCbE8dxYBiGv/goisKurTWordDh7U97MDkfEEFwOJywWq0oLS2Nui9e+vmAXC6DXC5Dfn4eGIbhx82QTh3y+WVlZRHP9dK9NvnDLj4WarVpvx8cy6GwQA6nh8bwMGAyAVotB4YRYXpag5aWHKjVLDweNxwOB0ymwKhZUo2WiTbKaPAzLD7uGsbgRED6t6woj38uak11bi7YoD1yPCRSAaZWq9P1BjcAGA/6/xSAskTemBEyO51OdHR0oKamJqXoXrK9q0ajEf39/UssZTyBgkRBjnP69Gn4fD60traG/KgBaSE//9pg5OWq8KWL2nBqcAJ/+ttx2B2OxYh16jc3sVgMtToXanXuYoGKHn6/HxRFYX5+nnfHo1mH4JRWaOrrjNUmEXO7yxtQK7HJkJPjXfz8wA3G6aSg1YqgUuVApTozatblCuTKaZqGWCxGbm4u5PL4c60iIdkbgnHBgef/+il2bq7C3sY6yKSSZdMRczgcQqSmZAhVfmQRcLfjIiNkHh8fR2NjY0oqhWS/negso7GxMRgMhojzq4SyzAzDYGZmBpWVlaivrw9ZG8MwfMlh9AATCwltx5VtWzCx4MPYrDHtNZHjzs/roVAokJeXB4oKeESBIJoJDONfDFypoFRGDlyFpr4C59Tn88PhcPBtmV6agUTGwutWQiT1gWEZsCywtHX3zKjZvLw8zM3NQSqVwmazgaa9kMlkUCoDe+3wls5oSEWUgOOAjoEJDE/rcfHueqhl1LKQWSBlzlkAXwz6fyWA44m8MSNkJuM5UgFJT8U7+QzDoLu7G1KpNKrQfjQdsGTgdDoxNDQEtVqNTZs28Y+T/XE8IhNxheLiYjRs2IBWikLv2Az+/nkf3N5IddeJIVIADQgEZkjTActy8HjccLmcGBiww+FQQqeToLZWBokk8kXp9zPQ6+dRUFAIpVLBp742bWLQ2emG1SqGmFKgtNSNgoLYmQeKAnJyVIv64mfmWs3Pn2npVKniTaNMXWHE7vTg5ffbUapVomVTQp5qTMS7Lt1utxAB1zcB/BtFUToAIgDnA7g9kTeueG12OAiZYwUQ3G43Tp06hYqKipiifulaZhLo2rhxIxyOM/ljQmTy40a72JxOJ7q6urBp06YQTe/6mnJUlxbi3c96MTAxF/G9sRAY5K5fjNZHDyqKRIFZxO3tufjoIzFIRLuhwYaWFsti6itnMXAVEGjQ6w3Q6Yp5L4ekvjQaYM8eDg6HDyKRD7oiKXx+P2gfE7VgJdRFDp1rxTBMyDTKaC2dAcuc9CkKwdC0HuPzJlwrVWL7xoqUjxNvzyyEXhnHcfMURd0L4OPFh37EcVz0RoAgZCw1lSriFY6QNFdDQwPfqB4N6eyZJyYmMDs7i9bWVrhcLtjtdgChga5YRDabzfw+PpLrpVLIcc3eZgxWzeOvJ0/DleDUCrfbDZPJtDjIPb4WmsMBHDsmRuA0BNZ6+rQW552ngEzmhs1mg9fr5XvJS0tLokrfyGQUCgoCx/D6Ant0TY4SdqcboDiw7Jn9eOC8RCci2UsH9pjRWzpZVpipGLTPj7eO96B/fA4Xt9ZDm5u8tnUsN1vI8bUcxz0F4Klk37dqLXM4OI7DxMQE5ubmEk5ziUQiXrkkUbAsi76+Pvj9fj7Q5fF4wLJsxIh1JExPT2NmZga7du2Kqwm1eUMJKnX5+Pvnfegdm4n5WrvdDpvNhtLSsqhucjicTgpiMRB8SsViwO0Wo7AwB7m5ObDb7bBaF5CbmwOj0bhIJhJEk0UkJPnubtoHTa4KfoaF20uHtHSSVFj8Cz28pTPQHGK1WuH1emE2mxeLZ5YKMSSCYIs5PmfC79/4GOft2ISWLVVJHS/ennml51GvCTKzLIvTp0+D47glkeRYSNbN9vl86OjoQGFhITZu3Mj/MOQ4JGIcSzd5aGgILpcLu3btSnidSrkMV57XhG3VZXjnxGnYXe6w4wZmeNG0F2Vl5XzrZCLIz19qHTkOKCwM7HetVivcbjcqKsr57+X3M4tlnRb4fD4oFIrFIJoy4me7vDQoikK+Jgc2pwccx8Jms/Hu95l+7fj14wCx2oGWzpmZaeTk5MLjCaTiRCKK32sHtmLxzwXLciHZA5+fwfsdA+ifmMPV5zclbKVjkXUlVTkJVr2bTcpAS0tLk1YfSSYA5nA40NnZiU2bNoWUf5K7usViwdDQEIqLixcjx6HrIAG5nJwcNDU1pXQONpYX42tX7cUHHf3oGp7irZrBYIRIRKGkpDTp/aNMBnz5y348/4IEHncgAn3jjX4oFIDRaALLMigtLQ1Zr0Qihkajhkaj5jutnE4XLBYzxGJJkNUOnfhoc7ohl0lhMZvhdntQWlrGrzdawUoi5ylQlBIYfeP3B4QYLBYLX88e2Gsr+AKYcEQioSZHiV1bq6HJSS5gFatcWYhqw3Swqi0zGQq3bds2vv80GSS6ZyZ56qamppD9LQl0yWQynHvuubBYLJidnUVfXx/UanXI0O3Ozk5UVlam3TUjl0pwSdt2bK0qwxsfd2JgZDQg5KfVphwIKi/ncOf/64PPRzoAOej1hsWSRF3M44Z3Wvl80evHAQrTM3PgWBZ1G6vh9Hh5ixVesBIgd6K92sE3GgnUajXUaiLE4FksoLFAJBLze+3gAGowmeUyCfbUb0TzlqqoQ+JTwUpLBgGrmMxEnTPViZJAYm72+Pg45ubmluSpw/fH4bW4NpsNRqORb2ovLy+PG5BLBsVaFRp0chRrtmFUb0vbjaOogJVmWRZzc/PIyVHFnQoSCbHqx8mNr7i4GA63BzKpBBKxOCS4R2INIlH6vdoBIQYlP5GTaJCTMbOkzJTjWIjFIjRv3oBzd2yKKOubLlZaZQRYhW62SCTC1NQU5HJ52uqcscjMsix6e3vBMMySctN4gS6KoqDVakHTNAwGA5qbm+F0OtHX1weaplFQUIDi4uJFa5r8ubDb7eju7sb2hkBv7azRijePd8Nsc8R/cwz4/Qzm5+eg1eYhNzf9C4/UjyuVKszPz0MqlUIsFmN+PrR+vDBPC7eXXjL0PpFe7TNET0SIQQqNRgqNRgOOY+HxBIbDa2QU9mwsxOZSDcCmlt2IdzPNWuYw0DSNiYkJqFSqlPedwYi2Z6ZpGqdOnUJRURFqamqiVnTFCnRNTEzAaDRi9+7dkEqlKCgowIYNG8AwDEwmE6anp9Hb2xvijifSF0xGszQ1NfEXR1lRHr5yxXk43jOCE6dHwKZgpWnaB71+HoWFhYLu7RiGxfz8PHJzc6HRBLYo4fXjer0BuTlK5Gk14Chx1PMaXD/OshyMRmPQdI9krbYINRWl+ELzFtiNc9DpdHC73ejv74fP5+N7kbVabUIR7XjCgistGQSsIjKTpm/SuidEiD/SnpkEusIF/RKt6CKpK47j0NLSsuRCEIvF0AUNRrPZbDAYDBgbG4NUKkVxcXHU1sypqSnMzs5i165dS3LIErEYe5s2o66yBG8d74bBakv4PJAurYCckXAuJrH0eXl5S6xScP14IIjmhdlqBxgfGA6QL0bIIxUHcRwHozGwp8/PLwBFYYnVDrjnkRVW1DkK7G2sw7bqMlAUha75aeTk5KCoqIi/4VosFn4QIZmcUVhYGDWVmEiTxbokc7JEJDJCzc3NcLlci1MX00e4m20wGDAwMBA10BWPyESAsKCgANXV1Qnt6cj+sq6ublEjzIje3l74fD4UFhaiuLgYarUaIyMjcDqdcVNaJQUa/OPl5+LT3lEc7xmJG+ALRKEtSXVpJQKfLzDuNhFLHwiiKfhqNY5jwfp80EeoHwewWG8u5+V7yDFiKaywLCCXSdFWvxGt9TUhwa3wRguxWIyioiIUFRWB485Mojx9+jT8fj+v+xUsoJ+ImN9Z7WYHC9ETGSGPxxPSO5wOCJlJX7Ver18yRTLRQhCXy4Wuri5s3Lgx6blXBEqlEhs2bMCGDRvg9/thMpkwOTkJg8EAuVweUvsd73udu30TNi9a6VmTNeLrbDY7HA572l1a4SDlpMXFxYtR7ORAUSKIZXJs2lgNmvbD5rDD5XLCZDKBYRgolcqYkz/CFVYoBEpkz2nYCKVcCo5l4VskMGlfjeXaB0+iDFYQ6e/vR05ODn/DWoaOqbSwYmQmQnt5eXlobm7mSSREcwQBiYz39PTwBSfJBLoIrFYrent7sX379sWmgfQhkUhQUFCAqakpbNq0CRqNZok7HiBL9Eq3Qm0ubrn0HHzeP44POwfhX7TSHAdYrRZ4vV6UlpYlVWQSD4G51QElE1kCKpmx4PJ4IaJEKCkuwIJdgdnZucUhcIBePw+O45bUj4djY3kxLmzeikJtgEjEwyJ/Al1gPt77irc/Ds9aOJ1O/qbr9XoxMjLCW+3g62XdutnxYLfb0dnZGXG6o5Cifn6/HzabDTqdbkmgK1YPcjBmZ2cxOTmJlpYWQZVSSLNIbW0tfw6IW+l2u2EwGNDT0wOGYVBYWIiioqIlFxAQsCy7t9WgtqJ4UQTBDKPRCIBDSUmJoOWFROwvVMkkPbAcC6vNCYvJBF1RIaSLe9a8vDwwDAu32w2bbQFeLw25XLbojqtQWpiHC1u2oro0tP5AJBKFuMajo6PIyckBRVFgGIav4hOLowfiCCiK4uvHtVot5uYCN5uZmRn09fUhNzeXF9B3uVyCSmJRFHUdgP8fADWAvwL4fziOi0mMZd8zEyH6nTt3RryTCUVmh8OBU6dOQSaTYePGjfzjie6POY7D8PAwHA4Hdu3aJdjFCwSCfadPn0ZDQ0PEXK9SqURVVRXv9hHLYLfbodVqUVxcjIKCghC3L1+dg+v37cLLb70Hp0OGnFxN2t1GwTgjSVSWVrowHH4/g7m5OeTn50MmVyBXJYfd5QUQKMHMzQ3Uj3NcoKuLYv2o0ohRWyQF57HD4ZDzZA0Gx3EYHR2F1+tFY2Mj726TbjfyBwDvjsciN+nkCw5uOhwOGI1G/Pu//zteeeUV7Nq1C42NjWhubhaifzofwDkAfAD+AuBmAM/GegMVJ3+WcqVCQCHyzNtJ3bLNZkNTU1PUFkev14uuri60xhFXi4XgQFdXVxfOP/98fg2JEJm45gqFAps3bxbUuun1ev5mlmyKiOM4WK1WGI1GmEwmyOVyPjouFov5sldNfiH+euI0RmfjasAlhIUFG5xOJ0pLSwSVf4oWRFPIpKAoKqTfWyqRoLW+Bq3bNkIqEcPrDQgXGo1GuFwu5OXloaioCAUFBRCJRBgaGoLP51siJkFAXHBCbgJC6vDvqdfr4XK5UFNTE/G73HXXXcjPz8f09DTuuusu7Ny5M5lTEfMCoyjqEQATHMc9Gut1y+Jmk3LH3NzcuEL06VjmeIGuRIjs9XrR2dmJsrIyVFZWprSOaJicnIRer+dz08mCoih+RvHmzZvhcrlgMBjQ1dUFu92O4uJiaDQaqFUKXPfF3Tg9Oo2/f94HD51c5xjBmb03jbKyUkFvaoExPvMoKloaRCPr1aiUcHq9aKipwPmNdchRnnmdXC5HRUUFKioqwLIsrFYrDAYDhoeH+eaQHTt2RF0zIatYLIZ0USSfxE+I1eY4jk+TxotmMwyDa6+9Fnv37k331ISAoqhyANcDuCzeazNOZpLX3bhxI8rK4qs9pEpm0lkFIOVAl8PhQHd3NzZv3pxSLXg0cByHwcFBeL3eiLnpVKFSqVBYWIiZmRk0NTXB7/djfHwcDoeDd8e/cvl5+Ht7P4am5pNcc0CcgWXZxb23IEsGcCYaHi/vXaDNwYHmFhTnxQ46kjEy+fn56O/vh9/vh1qt5lOABQUFKCoqilkgEmyNg9tdSTCWpumQcUfhSDU1dfvtt+M3v/nNybCHb+M47hRFUa0IuNb/wnHcYLxjZYzMRFxucHAQjY2NCUeBU7n70zSNjo4O6HS6iPnfRAJdpPKqsbFR0HxhcDdVLEuRCiwWC98gQtZcWlrKWyqj0Qjz8DCq1Ark15Xh1JgetD/+jTLQqWWAWByI7ApJZKJmEisaXpSnxoXNW1FTVpTwcTmOQ19fH8RiMbZv3w6KolBdXc3PaiYNMqR4pKioKKq4Q7DVBgIztufm5tDQ0BB1r02UOZPFr3/9a/z6179esqekKOpCAL8E8CWO404lcqyMkXl4eBhGo3GJuys0yAC6zZs3h0jzELdarVbj5MmTKCwshE6nW0x9hF6dk5OTmJ+fj1h5lQ5I2WgmXPb5+XmMj4+jubl5SZSdWCoSXXU6nTAYDBD53fhscBoGBw2VKgcy2VLhAZblMD8/D6VSEVK0IQTcbg9MJmPUApYcpQLnN9ZhR21FUjc9juNw+vRpyGQy1NXVhbVzSiIGrU6dCvCDZArUanXEz3Q6nejp6eHnokWy2uSGIXBq6kkAV3IcN57oGzIWAAuIzWlScik/+ugjPmgVC3q9nq9jDj6R4ftjv98Po9EIg8EAp9PJN0JoNBoMDQVGyTQ0NAga3HG5XHzZaFFR4hYmEZC9d6xAYjT4fD581jOANz/uhNXuhEKhQE5OQHggQOQ55Oaq+TproUC6mUpLS5ZkBiRiMVq31aC1fiNkSVapcRzHBys3bdqU1E2ApmmYTCYYjUZ+a0KCaBKJBE6nE52dnWhsbIxKVJZl8dFHH+HrX/86Ojs7QwxKEghZNEVRSgB2AGNBD7/Lcdz/E/MgmSKz3+9POZAVj8xEYtdoNGLnzp1JBbpYloXZbMb8/PyiBVKitrY24UaIRECKTKLpf6UKki5zuVzYsWNHWjcfD+3D++19ONEzBKfTBbfbBb+fgVqdi/z8fEHTT07nGeH/8AmXDRsrsLexDrmq5HP4LMuip6cHOTk5qK2tTWuNLMtiYWEhsDUxm/lBcA0NDTEJevLkSfzgBz/AK6+8gurq6lQ/XpCNzKol83nnnRc1pdDT0wORSIT6+vqUAl1utxudnZ2oqqqCSqWCwWCAyWTie3GLi4vjandFw/z8PMbGxrBz505Bi0xIy6ZEIsGWLVsE23tPzJvw+rEO9A+PLvYoM3A6XQACAbacnOg6YInAbncs6paVhpSUVpUU4sKWrdDlp1ZRx7Isuru7oVarQ+oIhIDL5UJ7eztKSkpgt9vh9XqRn5+PoqIi5Ofn89dcR0cHvve97+HFF19MuBQ3CtYvmY8fPx5R64umabS3t0eUEEq0oitWwQaZxGA0GsFxHIqKiqDT6RIKiJG2SJPJhMbGxnRmDS0BwzDo7OxEfn5+Qg0eycDhcKC94xRsnAID00a+NiDQwuiCy+UCTfugVCqgUuVEFdOPBJvNBocjND9doMnFvpat2FiekjsKIEDkrq4uaLXaqHnfVEG2R9u3b+e9KtJlZTQaYbFY0N7eDr1ej1dffRWvvPIKtkQY7J4k1i+ZT5w4scR9JiWgW7duXbIHTWSqBBCoPhsfH0dTU1Pcgg2apmE0GqHX6+HxePgOp0iCAxzHob+/HwzDLPEW0gUJolVWViaU2ksGVqsVfX19fAR/xmDBW592w2wLlWnmOA5utwculxNutwdSqXRxCmT0iZRW6wLcbjdKSkoC+t1yGc5r3IzGTRVpnR+WZdHZ2YmCggJUVVWlfJxIICW2DQ0NUbMvHMfhjTfewM9//vPFACKFZ555Jl03f3WTmdTBpoL29nZs27aNJ1yiga5YpZmjo6P8vOZk98ZEcMBgMMBms4WUVAJAV1cX1Go1amtrBa+FJjXsQua9gTOpuPBouJ9h8HHXMD7rG40ogsBxgM9Hw+kMWG0g1B0HAoUmNE1Dp9NBKpFg92JwS55mCybxUEhfspAgRK6vr48ppzQwMICvfe1reOaZZ9DY2Lg4Xzv6XK8EsX7JfOrUKWzatAk5OTkYHR2FyWRKOtBFQIpJyF5TgIkDfLWR0WjkL9rNmzcL6lrbbDb09PQI2qlFMDc3h8nJySXnNBjzZltCIgjEHXc6XfD5AvOdxWIxSktL0bCxAhfs3AJ1CsGtSJ9z6tQp6HQ6wdN8Ho8HHR0dcYk8OjqKL3/5y3j66afR0tIi5BLWL5m7u7tRUVGByclJSCQSbNu2LaVAF2mz1Ol0grtkDocDXV1d2LBhA3w+32KRhZgPoKUjzWMymTA4OIimpiaoVMlPXogFktbauXNnXA+FYVh82juC4z0jcdtSOQ4wGg3w+xlUFudjc6kalSVFfO14Ojc6hmHQ0dGB0tJSVFSkPl4mEgiRt23bFjOvPjExgZtvvhm//e1v0dbWJugasNrJzLJs0tMkCLq7u7GwsIANGzYsIWH4eJhoIHOeMpHnJaNnwvOPgV5fAwwGA3w+H4qKinglkUTd79nZWUxNTcW0mqmAbDXsdjsaGxuT8lCMVjve+rQbc6aFqMc2GAwo0OTi+kv3oq6yhC/QIJkCiqL485FMhZ3f78epU6dQXl4ueMwgUSJPT0/jS1/6Eh5//PGE6h9SwPoks81mw8mTJ1FdXb0k3J9ooItYth07dgjeME76m5uammKmnoILVRwOB/Lz81FcXByS2gjH2NgYLBZLSvv6WOA4DgMDA3yALpV9Pcuy+Lx/HB91DfEiCOTYCxYzztleiysvPCfqdyNdTgaDAR6Phz8feXl5Ud/j9/vR0dGByspKlJaWJr3mWPB6vWhvb8fWrVtjSiTPzc3hxhtvxH/+539i3759gq4hCOuPzPPz8xgeHuY7g8gPmOj+GDgjitfU1JRyrjgSSKGK1WpNmmwsy8JiscBgMMBisSA3N5d3PyUSCU82n88neCUaiRnI5fIlpY6pwGJz4q1PezBtMENEUdBKGHyhZRvqNiUezSWpHoPBAKvVitzcXL5emrjjZFRQVVVVyIQRIZAokfV6PW644QY89NBDuPjiiwVdQxhWN5k5jgNNJzZ/mOM4jIyMwGKxYOfOnZienoZUKkVFRUVIz2msqYukM8nj8WD79u2CVjARRU4AS/bvyYLjONjtdj6AJpVK4fP5oNVqsXXrVkGj4cH5aSHzsRzHoXt4Cub5KdTVVKW1jw0+HyaTia8rn5+fx6ZNm1LWW4sGUquwefPmmMogRqMRN9xwA37605/iiiuuEHQNEbA+yEyEAIIDXZOTk+A4jpdFjRfoYhgGXV1dyM3NTbo+Nx78fj+6urqQl5e3RHpIiGO3t7dDJpPxeXkSQIuknpEMfD4f3+QhdNCIWM0NGzYI7v7a7XacOnUKUqkUHMeFDBRI12MhRK6rq4uZ6rNYLLj++utx33334dprr03rMxOEIBfViqpzer1edHR0oKysLCTQJRYHlCQSIbLH4xFszlO0Y1dVVWVkz3bq1ClUV1fzbiSJig8PD8PtdvMXcqRBdYkcu6amRnDLRo69cePGVJsKooKmaZw+fRrbtm1DUVERGIYJaV9Uq9UoKipCYWFh0tHxRIm8sLCAL33pS7j77ruXi8iCIWOWGQj88NFgs9nQ1dUVcSjc7OwspqenUVdXF7FlMfgYZLCckHOegDMjYjJxbBJp37JlS1RXj1zIBoMBCwsL0Gg0/GSMWFsIUvwQ69iRYLcDY2MUWBaoruYQKbhLor+ZKGIhN/ZoZAt3x0kasKioKG76jvS719bWxsxs2O123Hjjjfj+97+Pm2++Oe3vlARWt5sNRCdzsKhfcJqC7I/9fj8MBgP0ej3cbjffixysTkm0tDKRiyXRcKGFCoAzteHJdFRxHIeFhQX+QlYoFLw7Hl7y2t3dnXShic0GvPKKGGRXJBIB+/czCOYUqUaLl8ZJBYTI8faxwSBpQKPRCK/XG7Xc1ufzob29PS6RnU4nbrrpJnzrW9/CV7/61bS/U5JY/WSOJOoXHOgKH7sZKWJNSin1ej3sdjvy8vIWa4XdKfXzxsPMzAymp6cFj4YDAaFBcgNKp6iEiA0YDAY+f6tQKDA2NpbSDejECQpdXSKQnYTRGLDO+/YFCkVIgUwmqtGItY8XWY6F8HJbMt9Lo9HwgwtibQncbjduvvlmfPnLX8Ztt92W6ldJB2uLzEQ+RyaTYevWrSlVdBFhQI/HA4qioNFooNPplsjOpgJyoyFFFUJGw4FA4cHs7OySm1i68Hq9GB0dxezsbIjFjqSxHQ0ff0xhcFAEYrisVkCn43DJJSy/lYnVoJ8qyJZASGtP5nvNz89jamoKKpUK5eXlUavyPB4Pbr31Vhw8eBC33367oAHOJLB2yOzxePgqnvAC+USJ7PP5+CJ7EixbWFiAXq+H2WyGSqWCTqfjc7fJILh+W+j0EKm8stlsGblJBFeMiUQi3kIRL4Y0hMSKBM/NAX/6kxhqdcDFtlopXHklA43Ggr6+voxsZRJtbEgFJEtQXV0NtVrNF6uQ+V5E3M/n8+GrX/0qLr30Utxxxx0ZJ/JVV12F8vJy/Pa3vw1/avVHsymKgtVqRXd3N+rr65fsh4I1i2NdbGS/Fjz9AQhMPSBut8PhgF6vx/j4OGQyGXQ63ZI9ZSQE3yTSUIqICCIyBwA7d+4U/GIhc6paWlr4G1hJSQlKSkpC5GcHBweRk5MTtU66tBS45hoWHR2BANi557LIzTWhr28gosZYunC5XDh16lRG3HZC5KqqKv5aCZ/vNT09jS9/+cuwWCxoaWnBN7/5zYwT+c0330RHR4fgGZdgZNQyT0xMYHh4GM3NzSF39mQquiyWgHVI5ocnIgNkT0kE3cIvSqI4UlNTI3iVEdlWECUMoa39yMgIHA5HQnXWwXXSRqORjwRHOidAYG8/OjqK5uZmwcUYia6W0JJKwJnyzw0bNsT8Pf1+P7797W+jqKgIeXl5+Oyzz/DGG29kjNBOpxMXX3wxbrnlFnR3d2fMMmeUzGNjYyguLg5xe5Mh8szMDKampuLWQccCiXrq9Xq+KEOn0/HFKvX19YJHZ0nBRmlpqeDtekQIgWXZlOusgxtC/H5/iHLp/Pw8Jicn0dzcLHhwkQTSMknkeHXcDMPge9/7HjZt2oT7779/WfbI3/72t3HZZZfB5XLh2LFja9PNrqysDFEbSXbOk9PpxO7du9PaZyoUCt7NIkUZ3d3dcDgcKC8vh0gk4vfrQoDEB2prawUvqiD6ZwqFIq29ffg5MRqNGB0d5edi19fXC763J0TORCCN9DpXVFTEJfIPfvADVFZWLhuRn376aVAUhZtuuglPPfVURj8ro5Y5WDoo0UAXsZhKpVKQxoBwTE1NYW5uDtu3b4fNZoNer4fD4UBBQQF0Ol3S1VbBIBdsJqx9puqsCSYmJmA0GlFZWQmTyQSr1cqneNJVLiX572CxfqFAep3jtUiyLIsf/vCHyM3NxcMPPyxoM0sstLa2wmq1QiKRYGEhIKV0/fXX43e/+13wy1a/m03InKjYHpnzVF5eLng9MRlc53a7lzRiEPldvV6PhYUFaLVaPuWV6I9OpktkotCEuO3l5eUZCaAER9vJ9yUpnnSVS0lqa+fOnYJHxIlFLi0tjXleWJbFPffcAwD4r//6r2UjcjieeuqptetmA4n3INvtdvT09CRdhpgIiCyrQqHgx3sGQyQS8S14RBZIr9djcHAQubm5fMormusZa7pEushknTW5wZGxp8EXOUVR0Gq10Gq1qKurCxlSx3FcSENINCwsLKC3tzejRC4pKYlL5P/9v/83aJrGk08+uWJEXg5k1DI//fTTqK2tRXNzc8w9GGkuyKRVKykpSVoEjtQDE/ldhULBp7xIcCid6RLxQFJymbjBkUAax3HYtm1b0pMg4imXEuXPVEbXxkMwkWN5cBzH4Wc/+xlmZ2dx6NAhweMAAmL1u9kvvfQSnn32WfT39+Oiiy7CgQMH0NbWFuLKkVxpY2Oj4GkQQoZNmzYJEowi6R2i90W8DVKwISRSrbNOBLFmMyWLSMqlSqUSc3NzaGlpEdxTYVkWp06dQnFxccxMAcdxeOihhzA0NISnn35aUOWWDGD1k5nA7XbjjTfewJEjR3Dq1Cns27cPV199Nf785z/jlltuwa5duwQnQyyx+3RBRNhpmgZFUbzbqdPpBHEnw/WshQTZcuTm5mYk/z05OYmRkRHIZDK+UCXYk0kHhMjxpHY5jsOjjz6K9vZ2PPvss4J7TBnA2iFzMLxeL15++WXcdddd0Ol0aGlpwfXXX4+9e/cKdtL1ej1GR0fTbmiIhEhRZZqmodfrodfreSG/kpKSlAQGjEYjhoeHBR9vE7z2goICwavdgEC3GdHilslkIQ0h6SqXEvH7wsLCuER+4okncOzYMbzwwgsZnUAqINYmmQHgJz/5CXbu3Ilrr70Wf/vb33D06FF8+OGH2LNnDw4ePIh9+/al/CNMTEzAYDBkZA+byHQJkreN1b4ZDZlS5gTOqFyWlJQIXsgCnLkJtbS0RFy7x+Phz0uyyqXEE8rPz48pmcxxHA4dOoS33noLR48eFbzrLYNYu2SOBL/fjw8++ACHDx/Ge++9h5aWFhw8eBAXXXRRQhaKiOLRNI3t27cL7ranMl0ivH0zPz+fz2WHr4/EDpqamgTf3xGZn0yMuAGSL/8kNdIkxx9LuZQQOS8vL6438dRTT+GVV17BK6+8IrhXk2GsLzIHg2EYfPTRRzhy5AjeffddNDQ04ODBg7j00ksj7klJHXROTo7gGmCAMNMliEKnXq+H1Wrl2zfz8/MxPj4Op9OZ9pjWSCAqG5lIbQGBLc3Y2BhaWlpS8oRiKZeKRCJ0d3dDo9HELZR55pln8Nxzz+HVV18VPA22DFi/ZA4Gy7I4ceIEDh8+jLfffht1dXXYv38/rrjiCqjVatjtdvT19WWk0ATIzHQJohyi1+sxMzMDiUTCR9yFtMqktDSe7lWqmJ+fx8TEhGB13OHKpR6PB1qtFtu2bYtpaQ8fPozf/e53eO211wQvFV0mnB1kDgbLsujo6MCRI0fwl7/8BYWFhRgdHcX//M//CD37B0Bm97DBddYlJSX8BZxM+2YskH7hdBQ8YmFubg5TU1Nobm4WfFvAcRy6u7shl8uhUChgMBiiKpe+/PLLeOKJJ/DnP/9Z8KwFAU3TuPPOO/HOO++A4zg8+OCDuOGGG4T8iLOPzMH48MMPcdttt+Giiy7Cp59+iqKiIhw8eBBXX321IFYoU9MlgNhR5UTbN2OBtBlmIkcNnBFczBSRe3p6oFKpQsakhgcW33vvPUgkErz++ut4/fXXM3LDIpibm8OxY8dw4403YmBgAHv27IHBYBAywJolc3V1NSorK/lqpiNHjuDVV1+FRqPB/v37ce2116K4uDipPXQmp0sAydVZR2vfjJV7JsUmmehOAgJtqbOzs3Gr+lIBKWZRKBRLRhMFg2EYPPLII3jmmWcgk8mwZ88e/OY3v1m2fHJRURGGh4eF9ATObjJHA2mfPHr0KF555RXIZDLs378fBw4cQGlpaUxik4IKlUqVkUAaUaHcuHFj0sEo0r5JSiiLioqg0+lCUjuk2CQTMj9AoOOMTJDMFJHlcnncc/+3v/0N999/P1577TUUFRWho6MDu3btEnQ90fDf//3f+P3vf493331XyMNmyRwPHMdhYmICR48excsvvwyWZXHttdfi4MGDqKysDLlgSB5Wp9MJPsgbELbOOjy1U1BQAKVSybu+QhfKAIHUmdFoRFNTU0aI3NvbC6lUGre89IMPPsC//uu/4rXXXhN8MEE8PPjgg3j++efx+uuvC53iy5I5GXAch9nZWRw9ehQvvfQS3G43rr76ahw4cABisRgnTpzAhRdeKLh8EJDZOmuWZTE6OorJyUnIZDLk5eUl3b4ZDxMTEzCbzWhqahJ820F00iQSSVwif/zxx7jrrrvw5z//OSOZi1j453/+ZzidTjz++OOZ8HqyZE4Her0eL730En7/+99jYGAAN954I77zne9gy5YtgrrXVqsVvb29GWnMB86kh4iEL2nfNJvNCbVvxsP4+Dg/+TJTRBaLxdi8eXPM837y5En84Ac/wJ/+9KeYVWCZwCeffIL77rsP77zzTqY+IkvmdDE9PY2rr74ajz32GAYHB3H06FHMzc3h8ssvx3XXXYf6+vq0LmAiep+JOmsgEIyamZmJGFVOpH0zHshw9kwUs5CgJUVRcW+gHR0d+N73voeXXnopJMK9XHjyySdxzz33hEzE+OUvfynkdMgsmdMFx3GwWCwhe1ir1YpXX30VL774IkZHR3HppZfi4MGDSbc5ZjJHDZwp/0w0GBXcvimRSPjIeLT65ZGRETidzoyUxpKMAcdxcbXMuru78e1vfxtHjhzBli1bBF3HKkKWzJmG3W7Ha6+9hqNHj6K/vx8XX3wxDhw4gNbW1pgXONHTykSdNRDIgS8sLKTs+rrdbj4yHt6+SWR8ibyS0BF9MkebZdm4RO7t7cU3v/lNPPfcc2hoaBB0HasMWTIvJ4J7sjs7O7Fv3z4cOHAA5557Lm8ZCREyVWdN0m6EaEIcP7x9UyQSQSqVZizYNTQ0BL/fH1fdZGBgAF/72tfwzDPPoLGxUdB1rEJkybxS8Hg8ePvtt3HkyBF89tlnOP/887F//3689tpruPnmm9HW1pYRizYwMACGYVLWy453/P7+fjgcDkil0qTbNxM5/tDQEHw+X9z1j46O4tZbb8VTTz2VkTLdVYi1QebHHnsMTz75JD/E+j/+4z9WajhXRkDTNN5++23ceeedUKlU2LVrF6677jpceOGFgu2VSR5WLBYLHm0nxw/fwybTvpnI8YeHh+H1etHQ0BBz/RMTE7j55pvx29/+Fm1tbel+tbWCtUHmmZkZlJeXw+12o6GhAa+++ip27NiR7mFXFR599FFwHIfvf//7eP/993H48GF88MEHfE/2P/zDP6QczSYNGaRWORNE7uvrg0gkinqjiNa+mej0zeHhYXg8nrhEnp6exk033YRf/epXOP/889P6XmsMa4PMBOPj47jyyivx0UcfCS4Qv9KINBGDYRh8+OGHOHr0KN59911s374dBw8exCWXXJJw0QHDMHxjfiaE74nFl0gkcfO8we8h7Zsmkwk5OTkxp2+OjIzA5XLFDabNzc3hxhtvxC9+8QtceOGFaX2vNYi1QebBwUFceuml0Ov1ePzxx/GNb3wj3UOuObAsi08//RRHjhzhe7IPHjyIyy+/PGozBJGTLS4uzkh5aTK10LGOQaZvRmrfHB0dhcPhwI4dO2IeX6/X44YbbsBDDz2Eiy++OJ2vtVaxush8++2347PPPgt57NChQ9i5cyeAQF706quvxuOPP44LLrgglbWuC5Ce7MOHD+ONN95AVVUV9u/fj6uuuorvwiGdVRUVFRmR+SFthkqlUlDXPbh90+v1QiKRxG36MBqNuOGGG/DTn/5UyCKMiHjhhRdw9913QywW41//9V/xrW99K6OflwRWF5kTwX333Ye8vDzcddddQh52zYI04R8+fBivv/46iouLcckll+D111/Hr3/964yMogmW2s1UNRXpBS8sLAwRFghv37RYLLj++utx33334dprr83IWgjsdjsaGhrwySefQCwWo7m5GV1dXYIP90sRa4PMH374Ifbu3QuHw4GLL74Y//7v/44vfvGLSR3jpZdewk9+8hPY7XZcfPHF+M1vfrOapxOkBI7j8MEHH+DWW29FVVUVcnNzsX//flxzzTVJ92RHAyEymRmdCUSq5Q6fgGEwGKBQKPDggw/iX/7lX3D99ddnZC3BOHLkCF5++WX84Q9/AADceuut2L9/P2655ZaMf3YCEITMGR+888ADD6Cqqgq7d+/GzTffnDSRgcAd/Pjx4xgcHMT4+Dief/554Re6wqAoCp988gmeffZZfPjhh3j88cfhdDpx66234pprrsGTTz6J2dlZxLn5RgVRudRqtRkj8sTERMSmDJlMhvLycjQ3N6O1tRU0TePee+/FzMwMjh07hpGRkYysJxiTk5Mhqi6VlZWYnZ3N+OcuJzI+s+O1115L+xjBe5udO3fCYDCkfczViB//+Mf8v+vq6nD33Xfjxz/+Md+TTYKH11xzTcSe7GggkyAKCwsz1nGUaJuk1+vF7373O9x999246aabhG7yjwqapkPWJRKJ1p13t6ZG4s3MzODFF1/EVVddtdJLWTZQFIXq6mr88Ic/xPvvv4/nn38eKpUK3/3ud3HJJZfgP//zPzEyMhLVYpOoeFFRUcaIPDk5CZPJFJfIbrcbt9xyC7761a/iq1/9KuRyOa688spl6YQqKyvD9PQ0//+pqamMZAlWEquqnDNWRPzkyZO49dZb8cADD+DGG29czmWtSnAcx/dkv/jii7Barbjqqqtw4MABvviDEFmn02VkigUQIAXp3opFZI/Hg1tvvRUHDx7E7bffvuxVgPPz89i1axfa29vBsizOP/98dHV1ZaTHPAWsjQCYEHj//ffx/e9/H7///e/5VFcWoTCZTHj55Zfx4osvYn5+HhdddBE++OAD/OIXv8hYo0KimmA0TeMrX/kKLrvsMtxxxx0rVs771FNP4ac//SkA4OGHH8Z11123IuuIgLOHzA0NDfjLX/4iyLCz48eP45//+Z/x8ssvZ8xarTQmJydx+eWXo7i4GDabje/JFrITanp6GvPz83GJ7PP58I1vfAN79+7Fj370o3VVly8gzg4yu91uqNXqkHLGiy66CL/5zW+SPtaPf/xjdHR0oL29He3t7euWzK+++ipcLhduvvnmkJ7sgYEBvid79+7dKRN7ZmYGc3NzcYns9/tx2223oaWlBf/rf/2vLJGj4+wgs5CwWq18nfOxY8fWLZmjweVy4S9/+QuOHj2K7u5uvif7nHPOSTiyOzs7y0sVxXoPwzD47ne/i7q6Otx///1ZIsdGlsyp4mwlczCCe7I///xznH/++bjuuutw/vnnR1VHSYbIP/jBD1BaWooHHnggS+T4yJI5GuLViWfJHAqapvHuu+/i6NGj+Pjjj3HOOefg4MGD+MIXvsD3ZCc6W4plWfzwhz9Ebm4uHn74YcHVStYpsmROFVkyR4ff7w/pyd61axdKSkpgt9vx0EMPxSXyPffcAwD4r//6ryyRE8faKOdcL3jhhRewceNG1NXV4Xe/+91KLydjkEgkuOiii/DEE0/g1KlT2LJlC5577jkcP34ct99+O/70pz/B5XIteR/LsvjJT34CmqazRF4hZLyccz3AbrfjRz/6UUjHDRlKt55BBApJc8ann36Kw4cP49/+7d+wefNmHDx4EJdddhlycnLws5/9DGazGYcOHcoSeYVwVrrZyWKVd9wsO1iWRXt7Ow4fPow333wTNE1jy5YtOHLkyLqrd14mCOJmZy1zAjgbOm6SgUgkwu7du7F792488MAD+POf/4yLLrooS+QVRtYfSgBnQ8dNqhCJRNi/f39GZkFHwksvvYTGxkbU1NTgtttuA8Mwy/K5awFZMieAs6HjZq3gbOhtTxVZMieAyy+/HG+++Sb0ej3m5ubw0Ucf4bLLLkvpWF6vF0888cRqKvJfU/jWt74FlUoFqVS6rnvbU0F2z5wASkpK8POf/xznnXceAOCRRx5JuXVu69ataGlpgd1uF3KJZx1Ib/tbb7210ktZNchGs5cZVqsVHR0d+NnPfpbJeb9rHmdZb3s2mr0Wsd4GAGQKv/71ryM+TnrbDx8+nO1tD0OWzFmsKXz3u98VrLd9vSFL5izWDNxuN9+TTZBqb/t6RJbMaxA0TePOO+/EO++8A47j8OCDD+KGG25Y6WVlHEqlEn6/f6WXsWqRTU2tQZjNZlx00UUYGBjAa6+9httuuw0+n2+ll5XFCiMbzV4HKCoqwvDwMD+rKos1h2wLZBbAf//3f6OpqSlL5Cyye+a1jAcffBDPP/88Xn/99ZVeSharAPHc7CxWKSiK+hWAHAD/xHHcUrWALM46ZN3sNQiKos4FsJXjuG+kS2SKokQURb1NUdQARVH9FEVdLtAys1hmZMm8NtEMoJWiqKGgP6lOKucAfI3juC0A/l8APxdqkVksL7JudhY8KIq6HcAejuNuW+m1ZJE8sgGwLEBR1I8B3A3AACDrZq9RZC1zFjwoiroewAMA6rnshbHmkCVzFiGgKGoKQDPHccaVXksWySEbADvLQVFULUVRpYv/Pg+AJ0vktYnsnjmLPABvUBQlBqAHcPPKLieLVJF1s7PIYp0g62ZnkcU6QZbMWWSxTpAlcxZZrBNkyZxFFusEWTJnkcU6QZbMWWSxTpAlcxZZrBNkyZxFFusEWTJnkcU6wf8fNWeIb/p5uigAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 训练数据集\n",
    "X_train = np.array([[3, 3, 3], [4, 3, 2], [2, 1, 2], [1, 1, 1], [-1, 0, 1], [2, -2, 1]])\n",
    "y_train = np.array([[1, 1, 1, 0, 0, 0]])\n",
    "# 构建实例，进行训练\n",
    "clf = LogisticRegression()\n",
    "clf.fit(X_train, y_train)\n",
    "clf.draw(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 习题6.3\n",
    "&emsp;&emsp;写出最大熵模型学习的DFP算法。（关于一般的DFP算法参见附录B）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**  \n",
    "**第1步：**  \n",
    "最大熵模型为：$$\n",
    "\\begin{array}{cl}\n",
    "{\\max } & {H(p)=-\\sum_{x, y} P(x) P(y | x) \\log P(y | x)} \\\\ \n",
    "{\\text {st.}} &\n",
    "{E_p(f_i)-E_{\\hat{p}}(f_i)=0, \\quad i=1,2, \\cdots,n} \\\\ \n",
    "& {\\sum_y P(y | x)=1}\n",
    "\\end{array}$$引入拉格朗日乘子，定义拉格朗日函数：$$\n",
    "L(P, w)=\\sum_{xy} P(x) P(y | x) \\log P(y | x)+w_0 \\left(1-\\sum_y P(y | x)\\right) \\\\\n",
    "+\\sum_{i=1} w_i\\left(\\sum_{xy} P(x, y) f_i(x, y)-\\sum_{xy} P(x, y) P(y | x) f_i(x, y)\\right)$$\n",
    "最优化原始问题为：$$\\min_{P \\in C} \\max_{w} L(P,w)$$对偶问题为：$$\\max_{w} \\min_{P \\in C} L(P,w)$$令$$\\Psi(w) = \\min_{P \\in C} L(P,w) = L(P_w, w)$$$\\Psi(w)$称为对偶函数，同时，其解记作$$P_w = \\mathop{\\arg \\min}_{P \\in C} L(P,w) = P_w(y|x)$$求$L(P,w)$对$P(y|x)$的偏导数，并令偏导数等于0，解得：$$P_w(y | x)=\\frac{1}{Z_w(x)} \\exp \\left(\\sum_{i=1}^n w_i f_i (x, y)\\right)$$其中：$$Z_w(x)=\\sum_y \\exp \\left(\\sum_{i=1}^n w_i f_i(x, y)\\right)$$则最大熵模型目标函数表示为$$\\varphi(w)=\\min_{w \\in R_n} \\Psi(w) = \\sum_{x} P(x) \\log \\sum_{y} \\exp \\left(\\sum_{i=1}^{n} w_{i} f_{i}(x, y)\\right)-\\sum_{x, y} P(x, y) \\sum_{i=1}^{n} w_{i} f_{i}(x, y)$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：**  \n",
    "DFP的$G_{k+1}$的迭代公式为：$$G_{k+1}=G_k+\\frac{\\delta_k \\delta_k^T}{\\delta_k^T y_k}-\\frac{G_k y_k y_k^T G_k}{y_k^T G_k y_k}$$  \n",
    "**最大熵模型的DFP算法：**   \n",
    "输入：目标函数$\\varphi(w)$，梯度$g(w) = \\nabla g(w)$，精度要求$\\varepsilon$；  \n",
    "输出：$\\varphi(w)$的极小值点$w^*$  \n",
    "(1)选定初始点$w^{(0)}$，取$G_0$为正定对称矩阵，置$k=0$  \n",
    "(2)计算$g_k=g(w^{(k)})$，若$\\|g_k\\| < \\varepsilon$，则停止计算，得近似解$w^*=w^{(k)}$，否则转(3)  \n",
    "(3)置$p_k=-G_kg_k$  \n",
    "(4)一维搜索：求$\\lambda_k$使得$$\\varphi\\left(w^{(k)}+\\lambda_k P_k\\right)=\\min _{\\lambda \\geqslant 0} \\varphi\\left(w^{(k)}+\\lambda P_{k}\\right)$$(5)置$w^{(k+1)}=w^{(k)}+\\lambda_k p_k$  \n",
    "(6)计算$g_{k+1}=g(w^{(k+1)})$，若$\\|g_{k+1}\\| < \\varepsilon$，则停止计算，得近似解$w^*=w^{(k+1)}$；否则，按照迭代式算出$G_{k+1}$  \n",
    "(7)置$k=k+1$，转(3)  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
